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AN ACCURACY STUDY OF FINITE DIFFERENCE METHODS 
IN STRUCTURAL ANALYSIst 
By Nancy Jane Cyrus* and Robert E. Fulton** 

NASA Langley Research Center 

SUMMARY 

An accuracy study is made of central finite difference methods for solving 
boundary value problems in structural analysis which are governed by equations 
with variable coefficients leading to odd order derivatives.. Two methods are 
studied through application to beam-columns with nonuniform inplane loads and 
noriuniform stiffness. Definitive expressions for the error in each method are 
obtained by using Taylor series to derive the differential equations which ■ 
exactly represent the finite difference approximations. The resulting differ- 
ential equations are accurately solved by a perturbation technique which yields 
the error directly. A "half station" method, which corresponds to making finite 
difference approximations before expanding derivatives of function products in 
the beam-column differential equations, was found clearly superior to a "whole 
station" method which corresponds to expanding such products first. 

'The material included herein was carried out by the first author in par- 
tial fulfillment of the requirements for a degree of Master of Science in 
Mathematics at Virginia Polytechnic Institute. 

*Mathematician, NASA Langley Research Center. 

**Aerospace Engineer, NASA Langley Research Center, Member, ASCE. 
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INTRODUCTION 


The differential equations governing the behavior of beams, plates, and 
shells are often solved by approximating the derivatives by finite differences 
and solving the resulting algebraic equations on a digital computer. In design 
analyses of complicated structures, such as civil engineering shell structures 
or aerospace vehicle structures, the number of simultaneous equations resulting 
from finite difference approximations can be sufficiently large to exceed the 
capacity of the computer or introduce round-off error. For such problems, the 
accuracy of the difference procedure can be a critical item in obtaining mean- 
ingful design results. In reference 1, for example, it was found that accurate 
answers for the stress in a shell could not be obtained by using certain finite 
difference approximations unless the mesh spacing was smaller than machine 
capacity permitted. 

The most popular difference approximations are the so-called central dif- 
ferences which are given in textbooks on numerical methods. There are alternate 
formulations of central differences which can be used when odd order derivatives 
occur in the differential equations. Such a situation results in structural 
problems, for example, when inplane loads are not uniform (a column loaded by 
its own weight or a shell of revolution) or where the stiffness of the struc- 
ture is nonuniform (a tapered beam or a variable thickness shell) . 

The purpose of this paper is to investigate the accuracy of two alternate 
forms of central finite difference approximations used in the solution of struc- 
tural problems. A new approach, for studying the accuracy of finite difference 
or finite element methods is presented and utilized. The study is confined to 
beam-column problems j however, the approach and conclusions are applicable to a 
wide class of plate and shell problems. 



SYMBOLS 


El(x) bending stiffness of beam 

f(x) nondimens ional tension in beam or string 

g(x) nondimensional stiffness of beam 

h finite difference spacing 

L(y) linear differential operator 

N(x) tension in beam or string 

p(x) nondimensional lateral load 

q(x) lateral load 

x axial coordinate of beam or string 

y deflection of beam or string 

Y deflection function in perturbation series (see eq. (l2)) 


STATEMENT OF THE PROBLEM 

Consider a general beam-column (fig. l) with nonuniform stiffness El and 
nonuniform inplane load N (taken positive in tension). The well-known differ- 
ential equation governing the lateral deflection y of the beam is 

(Ely 11 )" - ( Ny ' ) ' - q(x) =0 (l) 

where q(x) is the distributed lateral load and where primes indicate differ- 
entiation with respect to x. This equation can be solved by finite differences 
by dividing the beam into stations of equal spacing h. The quantities El and 
N are known, and finite difference equations are written in terms of the dis- 
placements at the ith station (i = 1,2,3 • • •)» 

In the present paper, two different finite difference approximations are 
considered. For convenience one formulation is called the "Half Station" method 



and the other the "Whole Station" method. For the term (Ny ' ) ' in equation- (,1) 
these two methods lead to the following finite difference expressions: 
i. Half Station Method 


(Ny * )[ 


|[-( N y , ) i .i/2 + ®y'h+i/2 

^ pi -1/2^1 _1 - (%-l/2 + Ni + i/2)yi + N 1 + i/2y i+ l 


(2) 


or 


2. Whole Station Method 


&y" + N’y'] i = %yi_i - 2 yi + y i+ i) + §-(-yi-i + 


JL. 

h 2 


pi - ^-)yf-i - 2 %yi + pi + ~^y 1+1 


( 3 ) 


Note that the half station method is the natural result of making the finite dif- 
ference approximation before expanding the derivatives while the whole station 
method results from making the approximation after the expansion. The latter 
type of approximation is widely used (see, for example, refs. 2 and 3 ) • Corre- 
sponding choices for the term (Ely") are: 

1. Half Station Method 


(Ely")!. = 


1_ 

h 2 


(Ely" ) i _ 1 - 2 (Ely" ) ± + (EIy") i+1 


1 


= 4 { (Ei)i-iyi-2 - 2 


h^ 


(El)i-l. + (El)i y irl 


(El)i_i + Mfil)i + (El) i+1 


n - 2 


(El)i 


+ (El) 


1 + 1 . 


y i+1 + (El) i+1 y i+2 


(4) 
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2. Whole Station Method 


j^l y iv + 2(El) 'y* * * + (ElV’/'J. = 


(EI) ± 


i h 


2(H) 


2h- 


y±-2 - k y±-i + - ^i+i + yi+2j 

-y±-2 + 2yi-i - syi+i + yi+2] 


tEiyr 1 

+ “ p- 1 ' 271 + y i+1 

= Jj- jjcEI.^ - h(El)jJy ._ 2 + ^(El)i 

+ 2h(El)’ + h 2 (El)^jy 1 _ 1 + ^(EDi 
- 2h 2 (El)|y 1 + ^(El)i - 2h(El)i 
+ h 2 (El)JJy i+1 + jcEl)! + h(El)[jy :L+2 


( 5 ) 


While the preceding two sets of finite difference approximations are both of 
order they clearly lead to different coefficients for the simultaneous equa- 

tions in terms of the displacements at the ith station. Of concern here are the 
relative magnitudes of the errors in these different approximations. 


ERROR ANALYSIS AND RESULTS 


The usual approach in a finite difference accuracy study is to carry out the 
numerical solution to a number of problems for which exact solutions can be 
obtained and to compare the resulting numerical answers at each station with the 
exact answers. Such a procedure has the liability that comparisons can only be 
made for each problem at specific stations and the calculations must be redone 
each time the mesh size changes. 
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The approach used in this paper is one which has not been reported pbe- • , 
viously in the literature. The finite difference approximations are expanded in 
Taylor series. This procedure results in differential equations which are exactly 
equivalent to the finite difference approximations. The resulting differential 
equations are then solved by a perturbation technique and yield analytical expres- 
sions for the largest error term. These expressions are independent of mesh 
spacing, are directly comparable, and give a clear indication of the relative 
accuracy of the difference approximations not just at discrete points but over 
the length of the beam. 

There are two terms in the beam-column equation which are approximated by 
finite differences: (l) the nonuniform tension effect and (2) the nonuniform 

stiffness effect. It is convenient to consider these two effects separately. 

Effect of Nonuniform Tension 

To study the effect of the inplane, load term in equation (l) let El = 0. 

The resulting equation describes the behavior of a laterally loaded string sup- 
ported at each end and subjected to nonuniform tension. For convenience, the 
variables are nondimensionalized so that the length of the string is 1 and tension 
is 1 at the left end. This leads to the following problem: 

-(f(x)y’ )' - p(x) = 0 (6) 

y(x D ) = 0 y(x 0 + 1) = 0 

where f(x) now represents the nondimensional tension in the string, p(x) is 
a nondimensional lateral load, and x Q is the coordinate o,f the left end of the 
string. Application of the two difference patterns, equations (2) and (3), to 
equation (6) yields: 
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1. Half Station Method 


" ^ £i -1/2^1 -1 " ( f i-l/2 + f 1+1/2 Vi + f i+l/2yi+l “ Pi = 0 (7) 

2. Whole Station Method 

( hf 1 A 

f i + — 2~/yi+l " p i " 0 (®) 

Expand the finite difference recursion formula equations (j) and (8) about 
the ith point using such Taylor series expansions as: 

y i± i = ± hjTi' + 1 | 7*" ± . . . 

f i±i = f i *“i' + ^ f i" 1 • • • 

For both the half station and whole station method this procedure leads to a 
differential equation of the form 

L 0 (yi) - Pi + h 2 L 1 (y 1 ) + h J4 l2(yi) + • • • 55 0/ (9) 

subject to the boundary conditions 

y ± = 0 at x — x Q 

yq = 0 at x = x Q + 1 

The symbols L 0 , Lq, and L 2 are linear differential operators given by 

Lo(yi) = -(%yi , )‘ ( 10 ) 
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and 


1. Half Station Method 


/f.y.iV f.ty.tM f.tlty.t 

Ufa ) = - + hZi — + h, + h — 

1 ' ^ 6 8 24 . 


12 


L?(y .) „ _/W i + w + f i ,,y i iv + f i ,,, yi ,M + + £lVV 

1 l 360 120 96 144 584 1920 j 


(11a) 


2. Whole Station Method 

Ll(yi ) 

LgCyi') 





(111) 


Equations (9) and (.10) together with either (lla) or (lib) are clearly dif- 
ferential equations which represent exactly the finite difference recursion 
formulas. As h goes to zero, equation (9) approaches equation (6). The solu- 
tion to equation (9) , satisfying the appropriate boundary conditions, gives an 
analytical representation of the numerical finite difference answers. Unfortu- 
nately a closed form solution to equation (9) does not appear feasible because 
it contains an infinite number of terms. For a practical problem, however, h 
is perhaps 0.1 or 0.01 or even smaller. This suggests that equation (9) can be 
solved by a perturbation method. with the perturbation parameter taken to be h 2 . 

Let the solution y i to equation (9) be taken in the form 
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y± = Y 0 + h 2 Yi + 


( 12 ) 



Substituting equation (12) into equation (9) leads to 

+ . . 

subject to 


L 0 ( Y o) - Pi + b 2 


Lo(Yl) 


+ Li(y c 


. = 0 


(13) 


^oCxo) + k 2 

Y 0 (x 0 + 1) + h 2 



. . = 0 

+ . . . = 0 


If each order of error term is solved in sequence, the following series of prob- 
lems result: 


(1) 

L 0 (Yo) - Pi = 0 

Y o( x o) = 0, 

G 

>T 

0 

+ 1) = 0 {ih) 

(2) 

L 0 (Yi) + Li(Y 0 ) = 0 

Yi(xo) = 0, 

Y l(xo 

+ 1) = 0 (15) 


(5) . . . ■ • • • • • 

Note that since equation (6) is linear Y 0 given by equations (l4) is in 
fact the exact solution. From the form of yi it is seen that Yi can be 
interpreted as the first order error term in the finite difference results. The 
magnitude of Yi is therefore a measure of the error in the finite difference 
results as compared to the exact answer to the problem. A comparison of the error 
terms Yq resulting from two different finite difference approximations indicates 
the relative accuracy of the two approximations when the node point spacing is the 
same. 


Using this method, the error functions Yq corresponding to the half sta- 
tion and whole station finite difference approximations have been obtained for a 
family of problems. These problems are a string having a lateral load which is 
distributed uniformly and a tension force f(x) which varies as follows: 


9 



(1) f(x) = i for 1 < n < 6 

subject to the boundary conditions 

y(i) = 0 
y( 2) = 0 

and 

(2) f(x) = 1 + x n for 2 = n = 6 

subject to the boundary conditions 

y(o) = 0 


y(i) = 0 


For the case where f(x) is Linear (corresponding to f(x) = 1, x, or 
1 + x) the results for the half station and whole station finite difference 
approximations are exactly the same. In fact for f(x) = 1, both difference 
answers are the exact answer. For all other cases, however, the two difference 
methods lead to different results. It is useful to compare the results for the 

case f(x) = — in detail as a typical example. 

■x? 

For f (x) = i and y(l) = y(2) = 0 

x3 


Y 0 = 


x5 , 31 k _ 16 
5 75 75 


1. Half Station Method 


*1 


jjl x k + x 2 + 86 

1125 6 150 1125 


(16) 


(17) 
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2. Whole Station Method 


v 187 1,^4 3 31 2 26 

1 450 3 30 225 


( 18 ) 


A plot of the two error terms Yp over the length of the string is given in 
figure 2(a). Solutions were also obtained for the error terms in deflection for 
all of the remaining load functions f(x) noted previously; an additional plot 
of results, for the case f(x) = 1 + x3, is shown in figure 2(b). Detailed plots 
of the remaining solutions are not shown because figure 2 serves to illustrate the 
character of the results; an overall measure of the relative errors in the two 
methods will be shown later for all the solutions obtained. 

While errors in the deflections of the string are important, errors in 
numerically obtained derivatives should also be considered for a thorough error 
analysis. Therefore, results were obtained by using the finite difference 
answers for approximate curvatures (second derivatives). The second difference 
operator was applied to the difference results followed by Taylor and perturba- 
tion series expansions to yield: 

7i " = h^ 73 -" 1 “ 271 + 7i+l) 

= Y 0 M + h 2 Yp" + ||(Y 0 iv + h 2 Yp iv +...)+... 

or 

y.” = Y 0 ” + h^Yp" + j + . . . (19) 

The h 2 error terms in the curvatures for the two methods and for the case 

f(x) = — are as follows: 
x3 
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1. Half Station Method 


Yq" + 



164 

375 


x2 - x + 21 


( 20 ) 


2. Whole Station Method 


Yq" + 


Y iv 
x o 

12 


374 


2*1 x 2 + 6x - a 


75 


25 


( 21 ) 


A plot of the error in the curvature for each of the two methods is also given in 
figure 2(a) for this case and in figure 2(h) for the case f(x) = 1 + x3. Again, 
results for the remaining load functions will he shown later in the form of an 
overall measure of the relative error. 

Numerical calculations were also carried out for the deflections and curva- 
tures for the problems cited to determine if the analytical errors adequately 
represented the numerical errors. The data are not included here; however, for 
h less than about 0.1 all analytical errors agree with calculated numerical 
errors to within 1 percent. 


Effect of Nonuniform Stiffness 

To study the effect of nonuniform stiffness on the numerical results for the 
behavior of a beam-column, the tension N is set equal to zero and the difference 
approximations given by equations (4) and ( 5 ) are compared. Results are obtained 
for a simply supported beam having a uniformly distributed load. Here again the 
variables have been nondimensionalized to make the length of the beam and the 
bending stiffness at the left end each equal to 1. This leads to the following 
problem: 

[g(x)y’lT = 1 (22) 
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y(x G ) = o 
y"(x D ) = o 


y( x o + i) = o 
y" ( x o + i) = o 


where g(x) now represents the stiffness of the beam and the distributed load 
is 1. 

From equations (4) and (5) the two difference equations resulting from equa- 
tion (22) are 
1. Half Station Method 


^7^i-l y i-2 - 2 (%-l + Si) y i-1 + (Si-1 + 4 Si + Si+D y i 

- 2( §i + g i+1 )y i+1 + g i+ iyi +2 J = 1 

2. Whole Station Method 



(23) 


(24) 


As before, expanding y^ and about the ith point leads to the differ- 

ential equation 

L o(yi) + h 2 L i(yi) + ^ L 2 (yi) + • • • • = i (25) 

where, now 

Lo = [g(x)iy i , [] (26) 
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and 


1. Half Station Method 

Vi n- 1 -\r V 


IHv.l't rt iVv ■*» 


/ \ s± , y i v 7 „ 1v Si ’yi'" Si 

11 6 2 12 1 1 3 12 


L 2(yi) = 


. g i ,y i Vii ^g. V vi + g W + J_ iv iv) 

' Sl y i 12 144^ ** 


8o 

tv H » 


20 


360 


, S i Vy i“' ^ S i Yiy i” 
+ ■ + 


60 


36O 


(27a). 


2. Whole Station Method 


, \ S^^ Si'yV s i "y i iv 

Li(yi) ■ + -V- + 




Si^l 


viii gi’yi vii gi”yi vi 


80 


20 


360 


(2Tb) 


If solutions to equation (25) , taking into account (26) and either (27a) or 
(27h), are again taken in the form ('12'), the series of simpler equations (l4) 
and (15) are again obtained (with p = l). However, since the beam equation is 
fourth order rather than second, a boundary condition on bending moment must also 
be considered. The moment is taken to be zero at the ends of the beam; this 
leads to 


Y 0 " = 0 at x = x 0 and x = x Q + 1 


(28) 


and 


Y-," + 


l iv 

;x> 

12 


= 0 at x = x r 


and 


x = x 0 + 1 


for the zeroth and first order error problems, respectively (see eq. (19))' 
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Results have been obtained for 


g(x) = x n n = 2, 3, k 

and 

1 < x < 2 

for both the half station and whole station methods of approximating the deriv- 
atives. The error terms for both deflections and curvatures are shown in fig- 
ure 3 for the case g(x) = x? corresponding to the case of a linearly tapered 
beam. An overall measure of the relative error in the half and whole station 
methods is given below for all three cases. The analytical error results for 
both deflection and curvature also agree with numerical error calculations within 
1 percent for h less than about 0.1. 

Relative Errors of the Half and Whole Station Methods 
While results such as those given in figures 2 and 3 are usually sufficient 
to identify which of the two methods is superior for a given problem, identifica- 
tion of the superior method for specific results is sometimes difficult (see, 
for example, the curvature errors of fig. 2(b)). Moreover, a quantitative meas- 
ure of the relative accuracy of the methods is desirable. Probably the fairest 
comparison of their overall merit can be made by examining the root -mean-square 
values of the errors for the whole structure ; that is: 

*1 « 

for the error in deflection and 

Yl" = 
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for the error in curvature, where the integration is over the (unit) length of 
the string or beam. Thus, to assess quantitatively the relative merits of the 
half station and whole station methods for the various problems solved, the 
ratios 

Y l,half 

Y l, whole 

and 

yll 

1 l,half 

Y" 

1, whole 

have been calculated for each problem. The results are shown in figure 4. 

DISCUSSION OF RESULTS 

The results given in figure ^(a) show that for all problems studied, the 
error in the deflection resulting from use of the half station method is less 
than the error due to the whole station method - in some cases, by an order of 
magnitude. The investigation of the accuracy of the curvature approximations 
gives the same result in general. Thus, the half station method is generally 
superior for calculation of both deflections and bending curvature for the prob 
lems studied. 

While the results are a clear victory for the half station method, one 
exception occurs: for the case of the string with the load f(x) = 1 + x 2, 
the error in the curvature is 25 percent greater with the half station method. 
Curiously, the difference between the two methods is seen to be generally less 
in calculating the second derivatives of deflections than in calculating the 
deflections themselves; moreover, differences in the comparative error from 
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problem to problem are noticeably less with the second derivatives than with 
the deflections. Both of these results are unexpected. 

It should be noted that the analytical representation of errors in the 

present paper shows clearly the danger of using numerical data at a single 

station or a few points to characterize the error in a problem. A typical case 

is shown in figure 2(a) for f(x) = — . If comparisons are made of the curva- 

x3 

ture near the end x = 1, the whole station method appears much more accurate 
than the half station method; whereas figure 4(b) shows clearly that the average 
error with the whole station method is over twice as great. 

It should be noted also that the present approach to error assessment may 
aiso be useful for comparison of different finite element structural approxi- 
mations. In fact, the recursion formulas given by the half station method 
(eqs. (2) and (4)) are the same recursion formulas which occur for a finite 
element model consisting of rigid bars connected by rotational springs, which 
often is used to replace the beam-column of figure 1 (see, for example, ref. 4). 
Thus, the results of the present paper verify that the finite element model 
of reference 4 is a good representation of beam-column behavior. 

Reasons for the superiority of the half station method are not altogether 
clear, but may include the symmetry of the matrix of coefficients in this method. 
By contrast, the matrix of coefficients associated with whole stations is not 
symmetric. Matrix symmetry can be of great value for many numerical procedures 
associated with eigenvalue routines and simultaneous equation solving routines 
and, in some cases, is required for an efficient numerical solution of a large 
order system. 
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CONCLUDING REMABKS 


A new procedure has "been developed to determine an analytical representa- 
tion of the error in a finite difference solution and to allow a direct compari- 
son between two difference methods which is independent of mesh size. This pro- 
cedure appears to have considerable merit for assessment of the relative accuracy 
of finite difference and finite element numerical techniques of structural 
analysis. 

Using this procedure, a comparison has been made of the accuracy of two 
different finite difference methods for solving structural problems through 
applications to a spectrum of beam and string problems having the characteristics 
of nonuniform stiffness and inplane load. The methods investigated were a "half 
station" method which corresponds to making the finite difference approximation 
before expanding the derivatives of function products and a "whole station" 
method which corresponds to expanding such products first ; both methods are in 
use. It was found that, for the same number of stations, the average error in 
calculated deflection resulting from use of half station difference approxima- 
tions was always less than the error which would result from the use of whole 
station difference approximations. In some cases this error is reduced by an 
order of magnitude. The investigation of the accuracy of the curvature approx- 
imations gave similar results in general. Thus, the half station method is 
indicated to be clearly superior to the whole station method and its use in 
finite difference solution of structural problems is recommended. 
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(a) f(x) = l/x3. 

Figure 2.- Finite difference error in deflection and curvature for a uniformly 
loaded string with nonuniform tension, f(x) . 










Figure J.- Finite difference error in deflection and curvature of a uniformly 
loaded simply supported beam with nonuniform stiffness, g(x) = x3. 




Figure 4. - Ratio of root-mean-square errors for half and whole station methods. 


